#!/usr/bin/env julia
"""
Final test of GSI diagnostic file generation
"""

using Printf
using Dates
using Random

# Load the unified Diagnostics module
include(joinpath(@__DIR__, "..", "src", "Diagnostics", "Diagnostics.jl"))
using .Diagnostics

println("="^80)
println("GSI Diagnostic File Generation Test")
println("="^80)

# Create output directory
output_dir = joinpath(@__DIR__, "..", "results", "diagnostics")
mkpath(output_dir)
println("Output directory: $output_dir")

# Create test observations
Random.seed!(42)
n_obs = 1000
observations = ObservationRecord[]

println("Generating $n_obs synthetic observations...")

for i in 1:n_obs
    obs = ObservationRecord()
    obs.lat = rand() * 180.0 - 90.0
    obs.lon = rand() * 360.0
    obs.pressure = rand() * 1000.0 + 10.0
    obs.obs_value = 1013.0 + randn() * 5.0
    obs.obs_error = 1.0
    obs.background_value = obs.obs_value + randn() * 2.0
    obs.analysis_value = obs.obs_value + randn() * 1.0
    obs.ob_minus_background = obs.obs_value - obs.background_value
    obs.ob_minus_analysis = obs.obs_value - obs.analysis_value
    obs.qc_flag = rand() < 0.95 ? Int32(0) : Int32(1)
    obs.use_flag = obs.qc_flag == 0 ? Int32(1) : Int32(0)
    obs.obs_type = Int32(120)
    obs.obs_subtype = Int32(181)
    obs.station_id = @sprintf("STN%05d", i)
    obs.time_offset = 0.0
    obs.inverse_obs_error = 1.0
    obs.variational_qc_weight = obs.use_flag > 0 ? 1.0 : 0.0
    push!(observations, obs)
end

println("  Created: $n_obs observations")
println("  Used: $(count(obs -> obs.use_flag > 0, observations))")
println("  Rejected: $(count(obs -> obs.qc_flag > 0, observations))")

# Create diagnostic output
analysis_time = DateTime(2018, 8, 12, 12, 0, 0)
grid_size = (95, 57, 32)

diag = DiagnosticOutput(output_dir, analysis_time, grid_size)
diag.observations["ps"] = observations

# Compute statistics
stats = Diagnostics.DiagnosticFileFormat.compute_innovation_statistics(observations)
diag.innovation_stats["ps"] = stats

println("\nInnovation Statistics:")
@printf("  Mean O-B: %.4f\n", stats["mean_innovation"])
@printf("  RMS O-B:  %.4f\n", stats["rms_innovation"])
@printf("  Bias:     %.4f\n", stats["bias"])

# QC stats
diag.qc_stats["ps"] = Dict(
    "total" => length(observations),
    "used" => stats["n_used"],
    "rejected" => stats["n_rejected"]
)

# Convergence info
diag.iterations = 10
diag.converged = true
diag.cost_history = [100.0, 75.0, 55.0, 45.0, 40.0, 37.0, 35.5, 34.8, 34.5, 34.3, 34.2]
diag.gradient_norms = [10.0, 7.5, 5.0, 3.5, 2.5, 2.0, 1.5, 1.0, 0.7, 0.5, 0.3]

# Write diagnostic files
println("\nWriting diagnostic files...")
Diagnostics.write_all_diagnostics(diag)

# List generated files
println("\nGenerated files:")
for file in sort(readdir(output_dir))
    filepath = joinpath(output_dir, file)
    @printf("  %-45s %10d bytes\n", file, filesize(filepath))
end

println("\n✓ Test complete!")
println("="^80)
